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A secure  digital  multiplexing  technique  is  presented  which  uses  a 
pseudo-random  process  with  a fundamental  cycle  length  N.  The  method 
can  handle  D < N signal  sequences  of  length  M < N.  The  SNR's  of  the 
de-multiplexed  signals  is  N/4DM  for  large  N,  even  if  the  additive  noise 


power  is  large.  The  de-multiplexer  must  compute  the  discrete  Fourier 


transform  of  the  pseudo-random  sequence  that  is  used  to  encode  the  signals 
Information  is  carried  in  the  phase  differences  between  the  frequency 


components  of  the  multiplexed  signal.  De-multiplexing  requires  at  least 


DMNlog-N  arithmetic  calculations 
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Quadratic  Multiplexing — A New  Method  for  Secure  Communication 


Melvin  J.  Hinich 

Virginia  Polytechnic  Institute  and  State  University 

Introduction 

This  paper  presents  an  information  coding  technique  that  increases 
the  channel  capacity  of  a communication  link.  This  technique  can  be 
used  with  any  reliable  system  for  transmitting  information  from  origin 
to  destination.  For  a computer  to  computer  transmission  link,  for  example, 
a number  of  messages  are  combined  into  a signal  which  is  transmitted 
by  standard  high  frequency  methods  to  a receiver  which  demodulates  the 
signal  into  a replica  of  the  original  multiplexed  signal.  This  signal 
is  then  unscrambled  to  yield  the  individual  messages. 

The  increase  in  channel  capacity  results  from  information  which  is 
carried  in  the  phase  differences  between  frequency  components  of  the 
broadband  multiplexed  signal,  and  thus  the  technique  is  inherently 
nonlinear. 

Information  is  encoded  using  a cyclic  pseudo-random  process. 

The  receiver  must  generate  the  same  pseudo-random  process  that  is  used 
for  encoding,  making  the  technique  secure  from  unauthorized  interception. 

I will  now  discuss  the  pseudo-random  process  before  going  into  details 

# 

about  signal  processing. 

1.  A Discrete  Pseudo-random  Process 

I will  limit  my  exposition  to  discrete  time  processes  and  digital 
signal  processing.  The  analogies  with  a continuous  time  system  are 
straightforward. 


mmmmmmmmmutsk 
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Let  {x(t  )}  denote  a discrete  pseudo-random  white  noise  process, 
n 

The  time  index  t - nf  \ where  n is  an  integer  and  f is  the  sampling 
n 8 s 

frequency.  Algorithms  for  generating  sequences  which  appear  to  be 
random  and  uncorrelated  are  well  known  by  now  (for  an  example  see 
Kronmal  [1]).  These  algorithms  produce  sequences  of  numbers  which 
satisfy  the  statistical  properties  of  a stationary  white  noise  process. 
One  of  the  most  widely  used  algorithms,  the  congruential  generator, 
produces  sequences  that  obey  a uniform  density  on  the  (0,1)  interval, 
although  it  has  a fundamental  period  which  is  determined  by  the  bit 
length  of  the  computer  registers,  and  the  "seed"  value  used  in  the 
algorithm. 

It  will  simplify  the  Fourier  transform  mathematics  if  {x(tn>}  is 

periodic.  Let  N denote  the  fundamental  cycle  length,  and  thus 

x(t  ) - x(t  ,,,)  for  all  n.  This  cycle  can  be  set  by  suitable  choice  of 
n n+« 

bit  length,  or  by  re-initializing  the  generator  at  n - N using  the 
same  seed  value. 

From  here  on,  the  x(t  ) values  will  be  treated  as  if  they  were 

n 

realizations  from  a finite  variance  stationary  white  noise  process, 

where  for  simplicity  the  expected  value  Ex(tQ)  ■ 0 for  all  n.^  That  is, 

for  all  n Ex(t  + t)x(t  ) - 0 for  t i 0 (t  an  integer  multiple  of  f *) 
n n s 

and  Ex2(t  ) * o2. 

Q X 

Define  the  discrete  Fourier  transform 
N-l 

X(f)  - V x(t  )exp(-i2xft  ) . (1) 

n n 

n-0 

The  principal  domain  of  X(f)  is  the  band  0 <_  f < fg.  Then  it  follows 
that  for  each  n " 0,1,..., N-l 
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X(tn)  * N tQ  X(fk)exp(i2i,fktQ).  (2) 

where  f,  * (k/N)f  are  equally  spaced  discrete  frequencies.  Moreover 
1C  s 

X(fk>  * X*(fN_k)  where  the  star  denotes  complex  conjugate,  and  X(f_k> 

- X*(f.).  Since  Ex(t  ) * 0,  EX(f.  ) - 0. 

K OK 

The  following  three  properties  can  easily  be  derived  from  the 

assumption  that  {x(t  ) } is  white.  For  each  k ■ 0,1,..., N-l 

n 

1)  ReX(fk>  and  ImX(fk)  are  uncorrelated, 

2)  The  variance  of  X(f^)  is  E { Xff k> | ^ * No^, 

3)  X(fk>  and  X(f^)  are  uncorrelated  for  k it. 

In  addition,  from  the  central  limit  theorem 

4)  {X(fk):k  * 0, . . . ,N-1 } have  a complex  Gaussian  joint  distribution 
in  the  limit  as  H + ®.  For  a discussion  of  the  complex  Gaussian 
distribution,  see  Brillinger  [2]. 

2.  Multiplexing  The  Signals 

Suppose  that  D digitized  signals  of  duration  Mfg*  are  to  be 

multiplexed  and  transmitted.  In  order  to  simplify  notation,  let  the 

time  unit  be  the  sampling  interval  f and  supress  the  unit  by  setting 

t - n for  each  n.  Denote  the  D as  follows:  {a. (0),...,  a,  (M-l) } , . . . 
n li 

{a^jCO) , . . . ,aD(M-l) } . Assuming  that  D and  M<N,  each  a^  signal  is 
convolved  with  the  lagged  product  sequence  {ox^x(n)x(n-d) }.  These  con- 
volutions are  summed  to  form  the  quadratic  multiplexed  signal  (Figure  1) 

A “f1  -2 

) a,(m)o  x(n-m)x(n-m-d) . 
d-1  m-0  a x 


y(n) 


(4) 
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The  signal  (y(n)}  is  used  to  modulate  a high  frequency  carrier.  This 
modulated  signal  is  then  transmitted  to  receivers  which  have  the  same 
pseudo-random  generator  that  is  used  to  multiplex  the  signals. 

Taking  the  Fourier  transform  of  (4),  we  have 
D i N-l  , 

Y(fk)  • V Ad(fk)N'i  E ax^X(fk_£)exp(-i2*f£d),  (5) 

d=l  £*0 

where  f k » ^ for  k - 0,  1,...,  N-l;  and 
M-l 

Ad(f  ) - E ad(“»)exp(-i2Ttfkm)  (6) 

m-0 

are  the  Fourier  transforms  of  the  signals.  If  prior  to  digitizing  all 
the  signals  were  bandlimited  at  f^,  there  is  no  aliasing  if  fg  = 2fQ 
and  the  quadratic  multiplexed  signal  will  then  have  the  same  bandwidth 
as  all  the  original  signals. 

For  certain  applications,  such  as  speech  transmission,  it  will  be 

better  to  set  M = N and  use  a much  higher  sampling  rate  than  fQ,  i.e. 

f »f  . At  the  receiver  end,  the  de-multiplexed  signals  must  be  smoothed 
s o 

by  a low  pass  filter  of  bandwidth  fQ  in  order  to  have  a high  output 
SNR. 

In  order  to  facilitate  the  statistical  analysis  of  the  de-multiplexing 

technique,  conceive  of  the  signals  as  realizations  from  Independent  (zero 

2 

mean)  white  Gaussian  noise  processes  with  equal  variances.  Let  o^  denote 
the  signal  variance  (power). 

3.  De-multiplexing  the  Received  Signal 

Suppose  that  the  receiver  digitizes  the  high  frequency  version  of 
the  multiplexed  signal,  yielding  a time  compressed  y signal  plus  noise 


* 
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denoted  by  p(n)  = y(n)  + e(n).  Once  again  the  time  unit  is  supressed,  but 
this  unit  is  much  smaller  than  the  original  sampling  interval  fg1  when  the 
signal  is  time  compressed. 

In  the  frequency  domain,  P(fk>  “ + where  (e(f^)}  is 

the  Fourier  transform  of  the  noise.  Even  when  N is  large,  the  P(fk> 
values  can  be  computed  in  a short  time  by  using  an  IC  module  programmed 
with  the  FFT  algorithm. 

The  receiver  must  have  either  the  pseudo-random  sequence  {x(n)}, 
or  its  Fourier  transform.  Suppose  that  the  receiver  generates  the  un- 
synchronlzed  sequence  {x(n  + x)},  where  x is  an  unknown  shift.  It  will 
now  be  shown  that  the  signals,  delayed  by  x units,  can  be  extracted  with 
a high  SNR  ratio  when  N/DM>>1. 

Denote  the  Fourier  transform  of  {x(n  + x)}  by  X^Cf^).  Clearly 
X^(f)  - X(f)  exp  i2itfx.  Next  define  the  weight  function  for  each  d 
and  k. 


W.(k)-(o  «)‘2  T X(f.)X(f.  .)exp(-i2irf  d). 

j«0  J K*J  -> 

As  is  shown  in  Theorem  2 (Appendix) , for  each  m « 0, . 

N"1  * 

s.(m)  - N 1 Y.  wj(k>p  (f.)exp(-i2xf  m) 
k-0 


N-l 


.,  M - 1 


- ad(m-x)  + ek. 


(7) 

(8) 

(9) 


-1  2 2 

where  the  are  approximately  uncorrelated  Gaussian  N(0,4N  (DMoa+0£)). 

In  other  words,  the  double  Fourier  transform  of  the  array 

{N~^o”2X(f  )X(f  )P*(f  )}  extracts  the  D signals  with  an  output  SNR  of 

X j X—J  K 

2 

N/4DM  when  a^<<N.  The  NxN  array  is  transformed  using  two  successive 


■ ■ 
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FFT's  to  yield  an  DxN  array  whose  dth  row  is  a^On).  When  M<N,  the  last 

N-M  columns  are  discarded.  The  first  row  (d«0)  is  2N  £ a,(m-T). 

d-1 

The  signals,  however,  are  shifted  by  t units.  If  the  receiver  can 
store  the  N complex  numbers  {X(f^)},  there  is  no  synchronization  problem 
since  expression  (9)  holds  with  t *■  0.  The  signals  are  extracted  with 
no  delay  if  the  receiver  can  store  or  compute  (X(f^)f.  Note  that  the 
transforms  in  (7)  and  (8)  can  be  efficiently  computed  using  the  FFT 
algorithm. 

For  applications  where  M * N and  fs>>f  , then  the  output  of  the 

Fourier  transform  in  (8)  must  be  shifted  back  to  the  (0,fs/2)  band  and 

then  smoothed  by  a low  pass  filter  in  order  to  increase  the  output  SNR. 

The  SNR  will  be  f /f  D if  the  filter  bandwidth  is  f . 

so  o 
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2 2 3 2 

Consider  the  following  example.  Let  o * 1,  o = 10  , M = 10  , 

3L  E 

and  D - 103.  The  SNR  for  each  signal  Is  -30db.  If  N = 5xl06,  then  the 

output  SNR  of  each  of  the  10  de-multiplexed  signals  is  lldb  (N/4DM  * 12% 

2 

and  a <<  N) . 
e 

Rather  than  computing  the  DM  array  (sd(m)}  by  using  DN  FFT's,  the 
double  Fourier  transform  of  X(f ^ )X(fk_^ )P*(fk)  can  be  computed  for 
d * 1, . . . ,D  and  m = 0,...,M-1  by  DMNlog2N  arithmetic  steps  using  the 
idea  behind  the  radix  2 FFT  algorithm.  The  computations  can  be  organized 
in  a sequential  manner  so  that  the  entire  NxN  array  of  the  triple  product 
does  not  have  to  be  stored. 

As  a final  note,  another  signal  can  be  transmitted  with  no  inter- 

M-l 

ference  with  the  D signals  by  adding  £ (n)x(n-m)  to  y(n).  This 

m=0 

linear  component  is  uncorrelated  with  y(n).  Thus  the  (n)  signal  can 
be  extracted  by  taking  the  Fourier  transform  of  the  sample  cross  spectrum 
between  {x(n)}  and  {p(n)L  The  SNR  is  the  same  as  that  for  the  other 

signals . 
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Appendix 


Theorem  1.  The  spectrum  S (f)  of  {y(n)}  is  ^ |A, (f)|  . The  spectrum 


d-1 


of  {p(n) } is  Sy(f)  + oe. 


Proof:  For  each  f,  there  exists  a sequence  k(N)  such  that 


f,  * -5-—  + f as  N-*».  As  is  shown  in  [2],  S (f)  = lim  N 1 1 Y (f  .) 


k(N)  N 


N-h» 


k (N) ' 


In  order  to  derive  this  limit,  note  that  EN  |X(f^)|^  - and  thus 


S (f)  * cr.  Moreover  by  properties  3 and  4 (Introduction), 


EX(fj )X(f^)X(f^)  * 0.  Applying  these  results  to  (5), 


EN 


'_1iY(fk)i2-  £ £ vfk>A*(fk>N"1*£1  ^ 

1C  “ ■ ' r K S K j.q  £„0 


r»l  s»l 


expi2n(f£s  - fjr)o^4EN'2X(fk_j)X(fj)X'(fk_£)X  (f£). 


(Al) 


For  large  N,  it  follows  from  properties  2,  3 and  4 that 


EN'2X(fk_j)X(fj)X*(fk_£)X*(f£)  - o*  if  j =*  l + k/2  or  j + t = k 


3o  if  j - £ - k/2 


(A2) 


■ 0 otherwise. 

Consequently,  most  of  the  cross  terms  in  the  second  double  sum  are  zero. 


Applying  (A2)  to  (Al),  N_1  \ Y(f  . ) | Z - £ |A.(f)|*  as  N~°  since 


k(N) 


N-l 


d”l 
k(N) 


Y expi2trf  (s-r)  * 0 if  r / s,  and  lim  N ^ Y expi2irf  (r+s)  = 0.  When 
j»0  J N-+~  j-0  J 


the  signals  are  realizations  from  independent  Gaussian  white  noise 


processes,  |A^(f)|  = Mo  . The  spectrum  of  (p(n)}  is  S (f)  + o2  since 

• * y e 


the  errors  are  independent  of  the  x’s. 


— 
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Theorem  2.  Expression  (9)  holds.  The  error  variance  is  approximately 
4N_1(DMo^  + o^)  for  large  N.  Thus  the  output  SNR  is  N/4DM. 


Proof:  First  consider  the  case  t = 0.  From  (5)  and  (A2), 

5 1 * D * 

EX(fj)X(f  )P  (fk)  = T Ad(fk)(expi2irfjd  + expi2irfk_jd) 

d=*l 

D * 

= Y.  A , (f,  ) (expi2irf  d + expi2irf,d  exp(-i27rf  .d) 
d«l  d k 3 * 3 


(A3) 


if  j j k/2.  When  j = k/2,  the  expected  value  is  3 > A, (f,  )expi2irf  . d 
N-l  d<L  d * 

Since  Y exp(-i.2irf  d)  = 0 for  each  d = 1,...,  N-l,  it  follows  fi 
j-0 

(7)  and  (A3)  that  for  t * 0 and  d>0. 


:rom 


* N-i  * 

EW,(k)P  (f  ) - N"  2 £x(f  )X(f  )P  (f  )exp(-i2uf  d) 

j“0  3 * \ 3 


A*(fk)  + 0(DM/N) 


(A4) 


since  N ^ A, (f  )expi2irf  . d is  of  the  order  0(DM/N). 


d-1 


‘d'V 


k/2 


Transforming  (A4)  into  the  time  domain  as  shown  in  (8), 


Es^On)  » ad(n>)  m = 0,...,M  - 1 (A5) 

D 

for  large  N.  When  d =*  0,  Es  (m)  » 2N  Y a,(m). 

0 d-1  d 

When  x i 0,  the  right  hand  side  of  (A4)  is  multiplied  by 
expi2xfjT  expi2irfk_^x  ■ expi2irfkx.  Thus  Esd(m)  « ad(m-i)  plus  an 
error  whose  distribution  will  now  be  derived  using  the  stochastic 
assumptions  for  the  signals. 
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The  next  step  is  to  derive  the  asymptotic  distribution  of  the 

errors  sd  - Esd-  The  term  B(j,k)  =»  N_1X(f  )X(fk_^ )p*(fk)  is  che  sample 

cross  bispectrum  at  between  (x(n) } and  the  received  signal.  As  is 

shown  by  Rosenblatt  and  Van  Ness  [3],  the  variance  of  B(j,k)  is 

NS  (f,)S  (f.  ,)S  (f,  ).  By  Theorem  1,  S (f)  - DMoJ  <■  a\  Thus 
x j x k-j  p k p at 

Var  B(j,k)  - No4(DMo2  + a2). 

x a £ 

For  large  N the  B(j,k)  variates  are  uncorrelated  in  the  triangular 
region  {j  - 0,...,[k/2];  k ■ 0,...,N-1>.  This  triangle  is  called  the 
principal  domain  of  B(j,k).  The  following  symmetries  relate  the 
B(j>k)'s  outside  this  domain  to  those  in  the  domain: 

B(k/2  + l,  k/2  - l)  * B(k/2  - i , k/2  + t)  for  k even, 

B((k-l)/2  + t,  (k+l)/2  - L)  - B ( (k+1 ) / 2 - t,  (k-l)/2  + l)  for  k odd,  and 

* 2 

B (N-j,N-k)  - B(j,k).  Thus  there  are  approximately  N /4  degrees  of 

freedom  for  the  B(j,k)'s  in  the  square  {j  * 0,...,N  -1;  k * 0,...,N  - 1}. 

Expression  (8)  can  be  rewritten  as  follows: 

N-l  N-l 

8.(m)«(No  ; > ^ B(j  ,k)exp[-i2ir(f  d+f,  m)  ] . By  the  central  limit 

X J-0  k-0  k 

theorem,  the  variance  of  s,(m)  is  4o~4N~2VarB( j ,k)  = 4N  i(DMo^  +0^).  The 
factor  4 results  from  the  symmetries  of  the  B(j,k)  in  the  double  sum. 
Moreover,  the  sd(“>)  are  uncorrelated  for  large  N by  the  orthogonality 
inherent  in  Fourier  transforms  (property  3). 
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Footnote 


1. 


Simple  whiteness  is  not  sufficient  to  obtain  the  asymptotic  results 
used  in  this  paper.  Assume,  in  addition,  that  Ex(nx)x(n2) x(n3)  = 0; 

Ex2(n],)x2(n2)  = if  n1  i n£;  and  Ex(nx)x(n2)x(n3)x(n4)  » 0 if 

nl  **  n2  **  n3  **  n4*  1,11686  conditions  hold  if  the  process  is  Gaussian 
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